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ABSTRACT 


A technique for improving the numerical predictions of turbulent flows 
with the effect of streamline curvature is developed. Separated flows, the flow 
in a curved duct, and swirling flows are examples of flow fields where 
streamline curvature plays a dominant role. A comprehensive literature review 
on the effect of streamline curvature was conducted in the present study. New 
algebraic formulations for the eddy viscosity \i t incorporating the k-e 
turbulence model are proposed to account for various effects of streamline 
curvature. The loci of flow reversal of the separated flows over various 
backward-facing steps are employed to test the capability of the proposed 
turbulence model in capturing the* effect of local curvature. The inclusion of 
the effect of longitudinal curvature in the proposed turbulence model is 
validated by predicting the distributions of the static pressure coefficients in an 
S-bend duct and in 180° turn-around ducts. The proposed turbulence model 
embedded with transverse curvature modification is substantiated by predicting 
the decay of the axial velocities in the confined swirling flows. The numerical 
predictions of different curvature effects by the proposed turbulence models are 
also reported. 
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CHAPTER I 


INTRODUCTION 


The tremendous improvement of computer capabilities in the past few 
years, including memory and speed, enables accurate numerical predictions of 
turbulent flows. Due to the closure problem of the governing equations for 
turbulent flows, numerous turbulence models have been proposed. The eddy- 
viscosity type of turbulence closure modeling has demonstrated a variety of 
good numerical predictions both qualitatively and quantitatively. Among them, 
the k-e model is the most widely employed isotropic two-equation model. This 
model has been extensively applied to different turbulent flow problems. 
However, the standard k-e model appears to be insufficient in predicting the 
complex turbulent shear layers, such as flows subjected to curvature and 
rotation. 

Flows with streamline curvature are particularly of interest in engineering 
due to their frequent presence in real life applications. There are several types 
of streamline curvature problems which may occur in the flow field. The 
separated flows would be considered to involve local curvature; flows in curved 
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ducts can be classified as a longitudinal curvature problem; and swirling flows 
typify the transverse curvature. The appearance of streamline curvature could 
change the structure of turbulent flow fields drastically. For example, 
turbulence intensity is enhanced with the application of concave curvature, 
while convex curvature inhibits the turbulent mixing. The change of turbulent 
structure will influence the mean flow field and vice versa. Consequently, any 
models do not include the effect streamline curvature will fail when streamline 
curvature occurs. 

Modifications to turbulence modeling are therefore necessary to account 
for the effects of various streamline curvature. The ad hoc change of modeling 
constants, however, to fit experimental measurements is not desirable, as it lacks 
physical rationale and generality. A large increase in computing time is not 
feasible for practical applications either. In the present investigation, new 
formulations for the eddy viscosity with the effects of different streamline 
curvature are proposed. Since the transport equations of the Reynolds stresses 
have a better description of the curvature effects, algebraic forms of the 
equations, approximated from the Reynolds-stress model, are essential. By 
extracting the extra strain rate and the main strain rate for the flow fields with 
different streamline curvature environments, new algebraic expressions for the 
eddy viscosity are derived from the algebraic Reynolds-stress model. In the 
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proposed eddy-viscosity formulations, the flux Richardson number R f plays an 
important role as it dictates the effects of various streamline curvature on the 
turbulence structure through the amplification or diminution of the eddy 
viscosity. Different coordinate systems are employed for various effects of 
streamline curvature in accordance with geometrical flow characteristics. 

The implemented eddy viscosity will be incorporated with the k-e model 
to predict the effects of streamline curvature on the turbulent flow field. A 
Navier-Stokes flow solver embedded with the standard k-e model is employed. 
A second-order finite differencing scheme for the temporal and spatial 
discretizations, incorporating a quasi-damping scheme which is deduced from 
second-order upwind difference concept, is adopted for the convection terms. 
The proposed models will be verified by comparing computationally predicted 
results with turbulent flows subjected to different streamline curvature effects. 
The examples for model verification are a) the flow over a backward-facing 
step, b) the flow in a curved duct, and c) swirling flows. The numerical 
computations are performed on the Cray X-MP supercomputer at NASA-Lewis 
Research Center. 
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CHAPTER II 


REVIEW OF LITERATURE 


During the last two decades, turbulence modeling has become the most 
productive, controversial, and rapidly improving subject for researchers in the 
computational fluid dynamics field. It started with simple boundary layer flows 
on a flat plate, and today numerous successes of the predictions on simple shear 
layer flows have been reported. However, the numerical computation of 
complex turbulent flows— such as three-dimensional flows, flows subjected to 
curvature and rotation, separated flows and shock boundary layer interaction— is 
still being investigated. 

The drastic change of flow characteristics due to the presence of the 
streamline curvature has been demonstrated by many researchers [l]-[4] based 
on classical laminar boundary layer theory study. The objective of the present 
research work is to explore the effects of streamline curvature and swirl on 
turbulent flows, which effects are much more complicated than their laminar 
counterparts. A review of literature is therefore conducted in two 

complementary ways— by 1) examining the experimental works to help 
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understand the physics of curved flows, and 2) studying and evaluating both the 
theoretical research and numerical techniques— to implement the turbulence 
modeling with the inclusion of curvature and swirl effect. The detailed 
information of the literature review in these two categories, and their 
summaries, are tabulated in Table A.l (for experimental work) and Table A. 2 
(for theoretical and numerical work). 

2.1 Experimental Investigations 

A variety of experimental investigations have been conducted in the last 
few decades. Bradshaw and his co-workers [5]-[10] conducted a series of 
experiments on convex surfaces, concave surfaces and curved ducts. They 
found that the turbulence intensity is reduced by the application of convex 
curvature which indicates a stabilizing effect, while longitudinal vortices are 
induced by the application of the concave curvature which has a destabilizing 
effect that enhances the turbulence mixing. Curved shear layers exhibited non- 
equilibrium behavior and required different formulation for the stabilizing and 
destabilizing effects. The rapid response of boundary layers to the presence of 
the convex curvature, as well as a fairly rapid recovery when the curvature is 
removed, has been experimentally demonstrated. Hence, Bradshaw et. al. 
suggested that the effect of convex curvature be modeled in Reynolds-stress 
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equations, at least partly as dependent on the rapid part of the pressure-strain 
term in the Reynolds-stress equations. It was further proposed that the apparent 
mixing length, increased due to the concave curvature, be modeled roughly 
proportional to the Richardson number Ri, which is 

Ri = 2(U « /R) (2.D 

3U e / d n 

where R is the radius of curvature, n is the normal distance away from the wall, 
and U 9 is the longitudinal velocity. It was noted by Smits, Young and 
Bradshaw [7] that the ratio of shear stress to kinetic energy is increased by the 
concave curvature. The increment of Reynolds stress is caused by the increase 
in the transverse production term in the Reynolds stress equation as the 
curvature occurs. Experimental data also revealed that flows recover more 
rapidly on the convex surface than on the concave surface as the curvature 
disappears. 

Gillis and Johnston [11] demonstrated the characteristics of flat-plate 
boundary layers recovering from a sustained convex longitudinal curvature. 
Since the radial pressure gradient acts to destroy the size of the largest eddies, 
the radius of curvature is proposed to be a scaling parameter. The experimental 
. results exhibited a slow recovery of Reynolds stresses from convex curvature 


6 


effect, which is contradictory to Bradshaw’s results. This is attributed to the 
stabilizing effect which permanently attenuates the turbulence length-scale and 
reduces the effects of the upstream condition. However, the near- wall layers 
are not influenced very much by the curvature. 

Streamline energy spectra for turbulent duct flow with small streamline 
curvature has been examined by Hunt and Joubert [12] in terms of a 
Townsend-type two-component turbulence model. Their result indicated that 
the flow was primarily affected by a direct change in turbulent shear stress 
through a conservative reorientation of the turbulence intensity components. 

Ramaprian and Shivaprasad [13] have performed extensive 
measurements of the surface curvature effects on turbulent energy balance and 
triple correlations. Turbulent energy production rate is reported to be 
significantly reduced by convex curvature and confined to a region very close 
to the wall; however, it is only slightly enhanced by concave curvature. 
Diffusion of both momentum and turbulent kinetic energy is found to be 
suppressed by convex curvature and strengthened by concave curvature. 

Effects of convex and concave surface curvature with artificially 
controlled constant pressure were tested by So and Mellor [14]-[15]. Reported 
test results are similar for curvature effects with and without adverse pressure 
gradients. Turbulent energy and its production rate, normalized by U t , appear 
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to be similar to the flat-plate data in the near-wall region; hence, no influence 
on the flow by surface curvature near the wall is concluded. A coherent 
structure of Gortler vortices is generated at the onset of concave curvature and 
then breaks up as flow moves downstream, due to the increase of turbulence 
level. 

Tani [16] performed an experimental test on concave-wall flow. He 
compared the data with Gortler’s small-disturbance theory because longitudinal 
vortices induced by concave curvature will intensify the growth of disturbances 
into turbulence in a boundary layer. Measured data reveal a spanwise variation 
having a definite wave number whether the boundary layer is laminar or 
turbulent, but the determination of the wave number is unavailable in Tani s 
paper. 

After inspection of individual surface curvature effects, curved duct 
flows, in which there is a combination of both curvature effects, are examined. 
First, an experimental study on the development of steady, laminar, 
incompressible flow in a curved pipe was conducted by Agrawal, Talbot and 
Gong [17]. The data indicate that at the vicinity of the curved pipe inlet, the 
inviscid axial velocity profile transits from a uniform distribution to a vortex- 
type distribution. 
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White [18] also investigated wall shear stress in a coiled pipe with 
various inlet Reynolds number. The results show that the friction force 
increases with the Dean number D n , defined as 



with r is the radius of the pipe. 

Strong secondary flows are discovered in S-shaped ducts with square 
cross section by Taylor, Whitelaw and Yianneskis [19]. In the first bend, 
secondary flows are larger in the laminar flow case than in the turbulent flow 
case due to larger inlet boundary layer thickness in laminar case. Taylor et al. 
reported that secondary flows reach their maximum values at the exit of the first 
bend, while the second bend produces the secondary flow effect in the opposite 
direction. 

Humphrey, Whitelaw and Yee [20] reported Laser-Doppler anemometer 
measurement of the mean velocities and Reynolds stresses in a 90°-bend duct 
with a square cross section. The locus of maximum velocity in laminar flow 
is shown to move more rapidly toward the outer wall than that in the turbulent 
flow. Low turbulence intensity and low level of anisotropy are exhibited at the 
wall with convex curvature, whereas high turbulence intensity and high level 
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of anisotropy are introduced at the wall with concave curvature. Some 
geometry with lower Reynolds number (laminar flow) was tested by 
Humphrey, Taylor and Whitelaw [29], and it appears that the secondary flow 
is already established at the entrance plane and persists downstream more than 
10 hydraulic diameters. 

Turbulent flows in a 180°-bend pipe and in a 45745° S-bend pipe were 
investigated by Row [21]. In the 180°-bend pipe, it is indicated that the 
secondary flows increase to a maximum and then decrease to a steady value. 
The flow in an S-bend pipe demonstrates that the secondary flows cause a 
complete interchange of fluids near the wall region and in the central core. 

Humphrey and Chang [22] explored the turbulent flow field in a 180°- 
bend square duct through L.D.V. measurements. According to the experimental 
result, they concluded that in the downstream straight section to a curved duct, 
turbulence diffusion and redistribution processes force the flow to erase all 
memory of the force imbalance acting on the flow in bend. However, a very 
weak secondary motion persists in the downstream straight section because of 
differences in the cross-stream gradients of the Reynolds stresses. 

Based on the present review, it is found that there are few experimental 
studies on two-dimensional turbulent flows in S-bend ducts, and the work by 
Butz [23] is one of them. Due to strong secondary flows in S-shaped ducts, 


10 


large aspect ratio is required to eliminate the lateral velocity component. Mean 
velocities and pressure recovery coefficients were measured in Butz s work. 

Brinich and Graham [24] investigated the turbulent flow and heat 
transfer in a two-dimensional 210°-tuming curved channel having an aspect 
ratio of 6, a radius ratio of 0.96, and adiabatic and heated walls; and operating 
at three different inlet velocities. The results show that at the start of the 
curved section an abrupt increase in pressure occurred for both the inner and 
outer walls. This is due to the force required to change the stream direction 
from rectilinear to curved flow. Another sudden pressure rise appears at 75 
station because of a change in the longitudinal vortex development. 

2.2 Theoretical and Numerical Methods 

Since the early parts of this century, researchers have realized the 
importance and complications associated with curved flows. Some researchers, 
such as Hawthorne [3] and Rowe [21], have tried to solve the problem through 
a quasi-inviscid theory. Meanwhile, laminar boundary layer theory incorporated 
with various techniques was employed more often in the curved flow field, for 
example, Tani [16] used small disturbance theory to examine the instability of 
the longitudinal vortices; White [18] proposed an empirical algebraic solution 
for the resistance coefficient due to the curvature effect; Agrawal, Talbot and 
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Gong [17] adopted an asymptotic expansion correction (developed by Singh 
[25]) to the governing equations of motion to compare with their experimental 
measurements; Stewart, Cebeci and Chang [26] obtained an approximate 
algebraic formulation by using a series expansion for three mean velocity 
components; and a similar expansion technique was employed by van Dyke 
[27]. Humphrey [28] and his co-workers [29] tried to use the finite 
difference method to predict the flow field. Although all the above methods 
have the advantage of simplicity and less computing time, none of them can 
capture the characteristics of curved flows very well, especially when the 
curvature is large. This is mainly because the mechanism of turbulence mixing 
does not exist in these proposed theories, and hence the amplification and 
attenuation of turbulence could not interact with mean flow field. In the 
remaining part of this chapter, we will concentrate on the turbulence modeling 
and its theoretical development. 

Bradshaw [30] performed a very extensive review of streamline 
curvature effects generated by surface curvature, swirling flows, rotating ducts, 
and/or spinning pipes. He also provided a method of distinguishing "simple 
shear layers from "complex" shear layers. A simple shear layer is defined as 
one where the simple shear, 3U/3y, is so much larger than any other rates of 
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strain that the direct effects of the latter on turbulence are negligible. Complex 
shear layers are recognized as perturbations of simple shear layers by the 
imposition of extra rates of strain or body forces, or by the interaction with 
other shear layers. Bradshaw [30]-[31] also deduced some characteristic 
parameters for flows with streamline curvature by drawing an analogy between 
meteorological flows, such as buoyancy effect, and curved-flow fields, such as 
radial pressure gradient. He proposed a correction to the apparent mixing 
length with small curvature effects based on the Monin-0 boukhov formula 

L = 1 - p Ri ( 2 - 3 ) 

'o 

where p is a positive empirical constant of order 10, and l Q is the length scale 
for zero curvature flows. 

A critical evaluation of various turbulence models performed by 
Nallasamy [32] reveals that the first generation turbulence modeling techniques 
developed based on simple shear-layer flows can no longer predict complex 
turbulent flows. Lakshminarayana [33] recently conducted a fairly thorough 
review of turbulence modeling techniques for complex turbulent shear layers, 
such as flows subjected to curvature and body rotation, separated flows, and 
vortex flows. The algebraic eddy viscosity and the two-equation models, with 
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constant value of C^, are concluded to be inadequate for the prediction of 
complex shear layers. Lakshminarayana recommended that a proper expression 
for C^, incorporated with a two-equation model, be used for two-dimensional 
flows with separation, curvature, or rotation. It is suggested that Reynolds 
stress models be employed for cases with very severe extra strains, large 
separation, curvature, or rotation effects, because of the large anisotropy. 

In the following review, papers will be summarized based on the 
sequence of zero-equation (algebraic stress) models, two-equation models, 
modified two-equation models, and Reynolds stress models. 

Hunt and Joubert [12] employed the length-scale model in equation 
(2.3), proposed by Bradshaw [30], to calculate a two-dimensional curved duct. 
No detailed comparison between experimental data and numerical prediction 
was reported. 

Towne and his co-workers [34]- [37] computed some curved-duct flows 
by using a parabolized Navier-Stokes (PNS) solver with a two-layer eddy- 
viscosity turbulence model embedded in it. In the outer region, the turbulence 
model of Cebeci and Smith [38] is used, while in the inner region either the 
model of Cebeci and Smith [38] or that of McDonald and Camarata [39] is 
adopted. In spite of less computing time required for the PNS solver, the 
streamwise marching technique has its own limitation (i.e., viscous flows with 
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no or small separation), while streamline curvature does include flow separation 
effect. Numerical predictions of turbulent cases are shown to be less accurate 
than those of laminar flows because of the failure of the two-layer turbulence 
model in recovering the effect of extra rate of strain. 

The Baldwin-Lomax [40] two-layer eddy-viscosity model was tested by 
Loeffler, Jr. [41] by applying it to an S-shaped diffuser. With the Beam- 
Warming [42] implicit scheme embedded, a large separation zone was 
predicted; but the numerical results were not validated by comparing them with 
experimental data. 

A simple mixing length model, in which the turbulent stress is 
proportional to the local strain rate of the mean flow, was used by Anderson 
[43] to compute duct flows with streamline curvature and streamline divergence. 
In the free stream, the length scale is selected to be the duct height, while van 
Driest’s [44] model is employed at the near-wall region. Detailed comparisons 
with measured data are not available in the Ref. [43]. 

Irwing and Smith [45] modified the Reynolds-stress model proposed by 
Launder, Reece and Rodi [46], with the assumption of local isotropy as well 
as local equilibrium, and obtained algebraic equations for the Reynolds stresses. 
The extra Reynolds stress production term is included to calculate the effects 
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of streamline curvature, according to Irwing and Smith’s observation that 
curvature effects on the Reynolds stresses are quite large. 

Gibson [47] introduced an explicit form for the length-scale function, 
where the influence of the wall on the fluctuating pressure field is modeled to 
account for the curvature effects. It was suggested in Ref. [47] that curvature 
effects can be included by the relatively small production terms appearing in 
each individual Reynolds-stress equation. 

An algebraic turbulence velocity scale for flows with curvature was 
developed by So [48]. The approximation was made based on the Reynolds- 
stress equations, in which the pressure-strain term is modeled by Mellor and 
Herring [49], with the assumption of local isotropy and local equilibrium. 
With the aid of two-dimensional boundary layer approximation, the velocity 
scale is derived to be a function of Richardson number Ri. 

Two-equation models probably are the most widely used models for 
complex turbulent flows. Chang, Han and Humphrey [50] applied the 
standard k-e model to a 90°-bend duct with square cross section. The 
discrepancies exhibited the failure of the standard k-e model to account for 
large-scale anisotropy in the flow. The standard k-e model in conjunction with 
a parabolized scheme was employed both by Patankar, Pratap and Spalding 
[51] and by Pratap and Spalding [52] to compute curved ducts and pipes. 
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Numerical predictions appear to underestimate the strength of secondary flows. 
Two suggestions were then made to modify the turbulence modeling for 
complex turbulent flows: 1) solving Reynolds stress equations, and 2) 
implementing the eddy viscosity based on approximate algebraic Reynolds- 
stress equations. 

Murthy and Lakshminarayana [53] compared the Baldwin-Lomax 
eddy-viscosity model with the standard k-e model, in a curved duct, 
incorporating in both models a space-marching, non-iterative algorithm. 
Although the standard k-e model is found to be superior to the eddy-viscosity 
model, neither model can produce quantitative agreement with the experimental 
measurements. 

A general algebraic expression for was deduced by Pourahmadi and 
Humphrey [54] to modify the k-e model to account for curvature effects. A 
new eddy viscosity was obtained by the combination of Bradshaw’s length-scale 
model in equation (2.3), and a functional C^, which can be written as 

F ( fa ) = c“ ♦ O.C, * a,C? * a, = 0 <2 ' 4) 

where a t , a 2 , a 3 = a,, a 2 , a 3 (PJe, T-Vr, 3Ui/3Xj). Santi [55] applied both this 
model and the standard k-e model to a 180° turn-around duct, with uniform and 
non-uniform inlet flow conditions, and compared the results. Based on various 
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computational results, Santi criticized Pourahmadi and Humphrey s model for 
an error in finding the root of the F(Vc^) equation (Eq. 2.4) and for additional 
difficulties which occur in specifying a selection criterion for regions with 
multiple roots. 

Launder, Priddin and Sharma [56] proposed that and C t2 be 
functional, and that an extra term be added to the energy dissipation rate (e) 
equation to include the effects of curvature. In this model, the modeling 
constant C* can be expressed as 


C; = C^exp 


-3.4 


( 1 + R t / 50 ) 2 


(2.5) 


where 1^ is the turbulent Reynolds number, and — 0.09 is the value of in 
the absence of the effects of streamline curvature. Another empirical coefficient 
C* 2 is modified to be proportional to a turbulent Richardson number Ri t (based 
on a time scale of the energy-containing eddies) and is defined as 


C t ‘ 2 = C t2 [ 1 - 0.3 exp { -R t 2 } ] ( 1 - C c Ri, ) (2 - 6) 

where C c is a constant with a value of 0.2, and C t2 is the modeling constant of 
the standard k-e model and equal to 1 .92. An additional term appears in the 
e-equation which is attributed to the curvature correction on the production term 
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of the e -equation instead of the decay part. This model was employed by 
Sharma [57] and provided reasonably good predictions. 

Rodi and Scheuerer [58] also compared this model with Gibson s 
algebraic stress model (Ref. [47]) in calculating the curved shear layers. The 
discrepancy between numerical and experimental results demonstrated the 
inadequacy of both models in capturing the curvature effects in general. 

An effort to include the anisotropy in the approximate algebraic Reynolds 
stress equations was conducted by Galmes and Lakshminarayana [59] in 
predicting three-dimensional shear flows over curved rotating bodies. In 
addition to using the same modified C* formulation as in equation (2.5), they 
implemented the production term of the e equation, based on their analysis, by 
modifying the empirical coefficient C*, to include the rotation effect, as 
follows: 

c;, - c„ /. « C„ [ 1 * 03 ( 1 - Ri c ) exp( -R? ) ] (2 ' 7) 

where C el is the modeling constant of the standard k-e model with a value of 
1.44, and Ri c is the Richardson number of rotation. A rather complicated model 
was proposed by these authors in which a modified k-e model, coupled with a 
set of six algebraic Reynolds-stress equations, needs to be solved 
simultaneously. In each of the six algebraic equations, the pressure-strain term 
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as originally proposed by Launder, Reece and Rodi [46] is noted to have been 
modified by Galmes and Lakshminarayana so as to include the anisotropy and 
near-wall effects caused by rotation and curvature. 

Warfield and Lakshminarayana [60] implemented the algebraic 
Reynolds stress model proposed by Rodi [61], to modify the Kolmogorov- 
Prandtl eddy-viscosity relation and produce an anisotropic turbulence model to 
account for the effect of rotation. Algebraic expressions for a vector form of 
Cjp which contains four components (C^,, C^, C^, and C^), were derived so 
that the Reynolds stress tensors are not necessarily aligned with the mean strain 
tensor. However, the formulations for the vector are considered to be 
relatively complicated. 

Naot and Rodi [62] obtained the algebraic equations for the vector 
modeling coefficient in the k-e model, composed of and C^y, to consider 
the anisotropic effect. The algebraic expression was derived from an 
approximated algebraic Reynolds-stress model by simplifying the Reynolds- 
stress equation proposed by Launder, Reece and Rodi [46], with the convection 
and diffusion terms neglected (i.e., in local equilibrium). 

A new eddy-viscosity model for swirling flow, representing a type of 
flow with streamline curvature, was developed by Kim and Chung [63]. With 
an assumption of weakly swirling flows, an expression for eddy viscosity was 
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derived through Rodi’s algebraic stress model [61]. A relatively good 
agreement was displayed by choosing the modeling constant (3 equal to 0.25. 
Nevertheless, according to Cheng [64], some inconsistency was found in Ref. 
[63], and it is not certain that the success of Kim and Chung s model is 
attributable to the inclusion of Richardson number or to the ad hoc change of 
the empirical coefficient C*. 

Wilcox and Chambers [65] demonstrated the prediction of the 
streamline curvature effects on turbulent boundary layers by using the e-co 
model. It was suggested that the streamline curvature primarily affects the 
equation of turbulent mixing energy e , while the equation of turbulent 
dissipation rate 0) remains unaffected by the curvature. A curvature correction 
term was then added to the ^-equation to account for the centrifugal effect. 

The Reynolds stress model proposed by Launder, Reece and Rodi [46], 
coupled with modeled energy dissipation (e) equation, was used by Gibson and 
Rodi [66] to predict a highly curved mixing layer. The Reynolds-stress model 
developed for plane flow is considered to have the advantage that neither 
modification to the basic closure hypothesis nor changes in the modeling 
constants are required to predict strong curvature effects. This method 
apparently not only is much more complicated and computing intensive but also 
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requires more modeling approximation for higher order terms. Some qualitative 
agreements with the measured data were reported. 

Hah and Lakshminarayana [67] investigated turbulent wake flows 
including curvature and rotation effects by using three different turbulence 
models: 1) standard k-e model, 2) e/algebraic-stress model, and 3) e/Reynolds- 
stress model. Due to poor representation and inadequacy of the production term 
of the e -equation for curved flows, different forms of this term in combination 
with those three models were tested as well. It seems that the k-e model with 
modified production term in the e equation can predict the streamline curvature 
effect very well but fails to capture the effect of rotation. It was suggested in 
Ref. [67] that the e/Reynolds-stress model or the e/algebraic Reynolds-stress 
model be employed to account for the rotation effect. 

Lilley [68] demonstrated the turbulent flow prediction results of several 
different turbulence models, namely mixing length model, k-k / model, Reynolds 
stress model, and algebraic stress model, in turbulent swirling jets. An 
analytical formula for the length scale was presented to account for swirling 
effect. An extra term was introduced in the k/-equation to include the effect of 
rotation, so that as the swirl number increases, k / will be enhanced. Fairly good 
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agreement with measured data was shown for the mixing length model and k-k / 
model. 

2.3 Conclusions 

Based on the literature survey, the following conclusions may be 
reached: 

1) Streamline curvature in the plane of the mean shear produces 
considerably large changes in higher-order quantities of the turbulence 
structure of shear layers: e.g., second order variables— Reynolds stresses, 
turbulent kinetic energy, etc. 

2) Turbulent mixing is inhibited by the presence of convex curvature; 
hence, the stabilizing effect will attenuate the Reynolds stresses and 
turbulent kinetic energy. 

3) Concave curvature has a destabilizing effect which will not only enhance 
the turbulence intensity and enlarge the length scale, but also induce the 
Gortler-type vortex structure, which the extra rates of strain become 
significant. 
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4) Streamline curvature effect could be caused by surface curvature, 
swirling flows, flow separations, or rotation of the whole system, and we 
consider a turbulent flow field with this effect a complex shear flow. 

5) The algebraic eddy-viscosity and the standard two-equation models are 
not adequate for the prediction of complex shear layers. 

6) Although it fails to predict the flows with large extra strains, Bradshaw s 
buoyancy analogy seems to be the simplest model to calculate the 
corrected length scale, through a correlation in which a Richardson 
number is involved to include the effect of streamline curvature. 

7) The k-e model, with the aid of modification of modeling constants 
(through the algebraic Reynolds-stress equations) or adding new terms 
in k- and/or e -equations, appears to be the most plausible approach 
regarding model complexity and the computing time. 

8) The k-e model, coupled with a set of algebraic Reynolds-stress 
equations, does have a better coverage of the physical characteristics of 
extra strains (from the effect of anisotropy); however, it is evident that 
more computing time and a more complicated modeling procedure are 
required. 

9) Even with a more complicated and time-consuming Reynolds-stress 
model, better predictions of curvature effects are not guaranteed. 
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10) In spite of the great number of turbulence models proposed for curved 
or rotating flows, most will show good predictions in only some 
particular geometries but will fail predicting others. 

11) For engineering applications, modifications to the turbulence modeling, 
allowing greater simplicity and less computing time, are feasible, 
desirable, and necessary. 
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CHAPTER III 


THEORETICAL APPROACH 


3.1 Governing Equations 

The equations of motion for a steady, incompressible, adiabatic flow of 
a Newtonian fluid are based on the following fundamental conservation laws, 
and will be expressed in tensor notation for universality [69]. The continuity 
equation (conservation of mass) is given by 



(3.1) 


The momentum equations, or so-called Navier-Stokes equations (conservation 
of momentum) are written as 


ou^i- + (3-2) 

H 1 ax. ax. ax, 

where Uj represents the three instantaneous velocity components in Xj- 
coordinate direction. Pis the instantaneous static pressure, p is the density, and 
the instantaneous stress tensor x S j is given by 
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*ii = M- 


3u : 




dXj ax i; 


(3.3) 


where |i denotes the coefficient of fluid viscosity. 

With the present computer capabilities, we know well the impossibility 
of solving the instantaneous Navier-Stokes equations for all the detailed 
fluctuating properties of a turbulent flow and with spatial resolution of finest 
eddies. Therefore, we can only hope to resolve the time-averaged quantities. 

^ A 

By using the Reynolds’ decomposition, Uj = Uj + Uj and P = P + p, and time- 
averaging procedure [70]-[71], we can obtain the system of governing 
equations for the time-averaged, steady, incompressible, adiabatic turbulent flow 
field, which can be written as 


3U. 


3X : 


i = 0 


(3.4) 


TT 3U i _ _ 3P + 3 
P j 35C 3Xj + 3Xj 


1* 


, \ 


3U; 3U } 

3x: + ax: 

V J 1 J 


-pu iUj 


(3.5) 


where U ; and P are time-mean parameters, and Uj and p are fluctuating 
components. We note that excluding the last term, pu^, in equation (3.5) gives 
the momentum equations for laminar flows. The term of -pujUj is found to play 
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the same role as the simple Newtonian viscous stresses in a laminar flow, 
hence, the fluctuating term is regarded as Reynolds stress or turbulent stress. 
For i = j, -pUjUj represents turbulent normal stresses, otherwise, the Reynolds 
shear stresses. Since the Reynolds-stress term appears in equation (3.5), 
additional equations are needed to solve the system of equations. A transport 
equation for the Reynolds stresses can be derived from equation (3.5) by 
multiplying it with u ; and taking a time average of the resulting equation [38], 
[72], i.e. 

C. = P . + Ps . + D n - e.. ( 3 - 6 ) 

u ij u >J '} 

where 


C ;j (convection) = U k 


P.. (production) = - 


Ps^ (pressure-strain) 


9u iUj 

ax k 


au. 


au. 


+ u : u L 


U;Uk ax ' Uj “ k ax„ 


p ax, ax, 
k J 1 J 
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D.. (diffusion) 
ij v ' 


du.u. 


dX L 


UiUjU|c v ax/ p v ' ,jk ' 


L_L + ±.(S jk u i + 5 ikUj ) 


(dissipation) = 2v 

and v is the coefficient of laminar kinematic viscosity, 8j k and 8 ik are the 
Kronecker delta functions. 

For i = j, by defining the turbulent kinetic energy as k = X UjUj/2, i.e. 
summing the normal Reynolds stresses, equation (3.6) becomes the transport 
equation of turbulent kinetic energy, which is 


^ 3uj 5uj ^ 


= ? +p k + D k - 


(3.7) 


where 


P r 

P k 


( convecton ) = U k 



3Uj 

( production ) = - au k — — - 

dX. 


. . p 5u i 

(pressure-strain) = — — 

P ax k 
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D. ( diffusion ) = - 


ax. 


+ _P 


- v 


ak 

ax„ 


e (dissipation) = v 


a Uj du. 
ax^ 


and the pressure-strain term P k is usually neglected on the evidence of measured 


turbulent energy balance [72]. 

It is obvious that every time we take time-averaging to introduce more 
equations, the more unknowns (higher-order turbulence quantities) appear in our 
system of equations. Turbulence closure modeling is therefore required to 
balance the numbers of unknowns and the governing equations. 


3.2 Turbulence Modeling 

As indicated in equations (3.4) and (3.5), the time-averaged Navier- 
Stokes formulation has more unknowns than the number of equations. The 
technique to model the new unknown term (pu s Uj) to close the system of 
equations is called turbulence modeling. For the last few decades, numerous 
turbulence models have been proposed by researchers for various types of 
flows. Detailed discussions and comparisons of various turbulence models are 
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described in the review papers by Nallasamy [32] and Lakshminarayana [33]. 
Despite a great variety of turbulence models which have been proposed, they 
can be classified into two categories: a) eddy -viscosity model , and b) Reynolds 
stress model. The eddy-viscosity models are constructed based on the 
Boussinesq’ s eddy-viscosity concept which probably is the widest-applied 
approximation for turbulence modeling, which draws a similarity between 
Reynolds stress and viscous stress, and is given by 


-pu iUj = p t 


9U; 91T 

ax" + ax 



(3.8) 


where }i t is named eddy viscosity or turbulent viscosity. Various methodologies 
have been employed to evaluate the eddy viscosity |i t . Based on the number 
of differential equations, in addition to the time-averaged Navier Stokes 
equations, used in each technique, the eddy-viscosity models are categorized as 
zero-equation (algebraic), one-equation, two-equation, and multi-scale models. 
As mentioned in the review of literature and suggested in Ref. [73], the 
standard k-e model, a version of two-equation models and proposed by 
Launder and Spalding [74], will be employed and implemented in the present 
study to predict the effect of streamline curvature because of its extensive 
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verification, wide application, and easy modification. A detailed formulation 

of the k-e model will be described later on. 

In the Reynolds stress models, the Reynolds stress tensors are computed 
directly by solving the coupled partial differential equations as shown in 
equation (3.6), in which higher-order terms appear and need to be modeled. It 
is obvious that the Reynolds stress models require not only more computing 
time and computer memory because of more differential equations and 
unknowns, but also more thorough understanding of turbulence physics to 
model these higher-order terms. However, the Reynolds stress models have the 
advantage of predicting complex turbulent flows such as flows subjected to 
separation, rotation and the effect of streamline curvature. Hence, in the present 
implementation approach, the essence of the Reynolds stress models will be 
adopted to correct the formulation of the eddy viscosity. 

With the aid of Boussinesq’s eddy viscosity formulation in equation 
(3.8), the mean momentum equation, i.e. equation (3.5), may be written as 

dU 4 _ _ l 3P + < 3 - 9 ) 

Uj axj *p 3X; + aXj _ p y dx i ax s J_ 

and p^ = p + p, is the effective viscosity. By using dimensional analysis, the 
eddy viscosity p, can be expressed as 


32 


H,~ P\l 


(3.10) 


where V t is the turbulence velocity scale, and / is the turbulence length scale. 
It was postulated by Prandtl and Kolmogorov and later adopted in the standard 
k-e model that 

l-*L . v.-yiT 0.11) 

e 

which will lead us to obtain 


m = PQ — 


(3.12) 


where is an empirical coefficient. In equation (3.12), two more unknowns 
are introduced and will require the solution of two partial differential equations 
for the turbulent kinetic energy and turbulent dissipation rate. This is why the 


k-e model is identified as a two-equation eddy-viscosity model. 

In the standard k-e model, the transport equation for turbulent kinetic 
energy as shown in equation (3.7) was modeled as 


c k = 


P r + D k - 


(3.13) 


where 
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c k = u k 


9k 

ax. 


( 




au : au. 


ax. ax ; 


au ; 


ax. 


p ax k 


p t ak 


a k ax k 
V k k j 


and as for the turbulent dissipation rate e, it is expressed as the dependent 
variable of a differential conservation equation which is highly empirically 


modeled and is given by 


U. 


9e 

die 


i a 

P a><k 


( \ 
9e 

o“ ax" 


V 


J 


r e ^ 
11 k p 


( 


V 


au. au k 
a5T + ax" 


au ; 

ax. 


(3.14) 


-C 2 — 


which physically represents the convection, diffusion, production and dissipation 
of the turbulent dissipation rate e, respectively. In equations (3.12)-(3.14), C^, 
c k , G t , C tl and C t2 are modeling constants, and some typical values of these 
constants in the standard k-e model are recommended by Launder and Spalding 


[75] in Table 3.1: 
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Table 3.1 The values of modeling constants in the standard k-e model 


C, 


c„ 

c, 2 


0.09 

1.0 

1.44 

1.92 

1.3 


According to equations (3.12)-(3.14), the standard k-e model apparently 
does not have the capability to account for the effects of streamline curvature. 
The modeling constant in the eddy viscosity formulation, as shown in 
equation (3.12), is empirically tuned for the simple shear layer. Meanwhile, 
there is no mechanism in the model which can either amplify the turbulent 
intensity or eddy viscosity in the presence of concave curvature, or inhibit 
turbulent mixing with the application of convex curvature. Therefore, the 
expression for eddy viscosity in the standard k-e model is considered to be 
inadequate to account for the streamline curvature effect. It is evident that 
modifications to the standard k-e model are necessary to include the curvature 
effects. However, the ad hoc changes in modeling constants are not desired due 
to lack of physical explication. The implementation to the formulation of the 
Reynolds stress should be reasonable and feasible. 
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3.3 Implementation Approach 

It was suggested by Muck, Hoffmann and Bradshaw [5] that the effects 
of curvature be modeled in the Reynolds-stress equations. Lakshminarayana 
[33] also indicated that the Reynolds-stress equations can interpret the physical 
phenomena very well. As a result, the implementation approach will start with 
the Reynolds-stress equations. Unfortunately, more higher-order unknown 
turbulent parameters are exhibited in the Reynolds-stress equations. Launder, 
Reece and Rodi’s [46] Reynolds-stress model is therefore adopted because of 
its well validation and wide application, which is given by 

G. = P . + Ps.. + D. - e.. ( 3 - 15 ) 

ij IJ ij ij ■ J 

where 


C. 


ij 



3u, Uj 


p «- " 


au au, 

u.u. L + u ; u,. 


rk 


ax. 


j k 


ax. 



u.u 


- 5 ii 


V 




Ik 

3 


- C, 


) 



) 
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d , = c _Li 
,J ‘ ax. e 


— au u — 
u U i — + u.u 

. m ax J m 


au kUj 

ax„ 


a UjUj 


+ ILU 

k m ax 



e 


and C, is an empirical constant with a value of 0.11. The modeling constants 
Cj and C 2 are inertial and forced retum-to-isotropy constants respectively, where 
their values will be discussed later on. It should be noted that in this model, 
a local isotropy is assumed, which will be valid as long as the cross flow 
(lateral or circumferential velocity) is not very large. Even the system of 
equations now is closed, solving these six partial differential equations coupled 
with k- and e- equations is still not feasible for today’s computer and 
engineering applications. 

An approximation of the Reynolds-stress model proposed by Rodi [61] 
is very physics-explicable and economical— the net transport of Reynolds 
stresses UjUj is proportional to the net transport of turbulent kinetic energy k 
with a factor of uju~/k, i.e. 


C ,- D , = 


(3.16) 
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In this approximation, it is assumed that UjUj/k varies but slowly across the flow 
field. By combining equations (3.13) and (3.15) with equation (3.16), an 
algebraic expression for uju” is obtained as 


Vi 

k 


1 

e 


P..-±5. P 

•J o IJ r 



where 


(3.17) 


This approximation will be invalid only when u ; Uj/k change greatly, i.e. D( 
ujuj'/k ) / Dt is not negligible, so it is not suitable for the near-wall region. 
Although equation (3.17) is in algebraic form, it is difficult to incorporate with 
the k-e model because the six Reynolds stress components are coupled. In 
order to further simplify this algebraic formulation, the understanding and 
analysis of flow characteristics are necessary. 

There are three kinds of basic streamline curvature: 1) local curvature 
such as separation bubbles and flows over an airfoil, 2) longitudinal curvature 
such as flows through a curved duct, and 3) lateral curvature such as swirling 
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flows. The coordinate system for each flow field is somewhat different, and 
also both the main strain rate and the extra strain rate vary in each case. 


3.3.1 Local curvature 

For flows subjected to local curvature, the streamwise direction basically 
does not change much, for example the flow over an airfoil or the separated 
flow in a backward-facing step geometry. In this case, the two-dimensional 
Cartesian coordinate system is selected. The streamwise direction is defined to 
be along the x-coordinate, i.e. i, j, or k = 1, and the transverse direction will be 
in y-coordinate as i, j, or k = 2. The velocity components are defined as U — 
U lf u = u t in the streamwise direction, and V = U 2 , v = u 2 in the transverse 
direction respectively. According to Bradshaw [30], the extra rate of strain is 
3V/0x besides the main strain rate d\J/dy for the flows with streamline 
curvature in Cartesian coordinates. By keeping terms associated with these two 
strain rates, and through some elaborate manipulation, the primary Reynolds 
shear stress -uv is obtained as 


-uv 


k 2 9U 2 
e 0y 3 



R 2 + 4R f + 1 
1 - R r ~ 


where the flux Richardson number R f is 


(3.18) 
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R dV/3x 

f au/ay 

The detailed derivation procedure will be described in Appendix B. By 
collaborating with the Boussinesq’s eddy-viscosity concept as shown in equation 
(3.8), an algebraic expression for the eddy viscosity p., can be derived as 


k 2 26 

H, = p T - 


P R f z + 4R f + 1 

1 - R f - <j) _L _ 1 

f e 1 - R, 


(3.19) 


With the flux Richardson number embedded in the eddy-viscosity formulation, 
streamline curvature can influence the eddy viscosity and so the Reynolds stress 
directly. It is evident that the flux Richardson number R f is positive for convex 
curvature ( stabilizing effect), and so the eddy viscosity is reduced. On the other 
hand, concave curvature ( destabilizing effect) gives a negative R f and will 
enlarge the eddy viscosity. This modification is therefore qualitatively 
consistent with the physics of the problem and conclusions reached from the 
experimental studies. The effects of streamline curvature vanish as R f = 0, and 
the eddy viscosity becomes 
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(3.20) 


k 2 2 , 


P 

e 


In order to determine two retum-to-isotropy modeling constants, C t and 
C 2 , in equation (3.17), the above equation will be matched with the correlation 
for the eddy viscosity in the standard k-e model at the free stream condition 
where there is no curvature effect and the production rate is equal to local 
dissipation rate, i.e. where local equilibrium prevails, and hence 

where 


The values of C, and C 2 have been determined empirically by several 
researchers such as Launder, Reece and Rodi [46] (C, = 1.5, C 2 = 0.4), Gibson 
and Launder [76] (C, = 1.8, C 2 = 0.6), Gibson and Younis [77] (C, = 3, C 2 = 
0.3), and Gibson and Launder [78] (C, = 2.2, C 2 = 0.55). Kim and Chung [63] 
also observed that 1.5 < C ( < 1.8, and 0.5 ^ C 2 ^ 0.8. However, the variations 
of Q and C 2 are expected not to affect the results substantially, which also will 
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be demonstrated in the numerical predictions later on. In the present model, C, 
= 1.5 is selected which will yield C 2 = 0.76 from equation (3.21). 

3.3.2 Longitudinal curvature 

The streamwise direction can change drastically for the flows subjected 
to longitudinal curvature, for example the flow in an S-bend duct or any curved 
ducts. Consequently, the two-dimensional Cartesian coordinates may not be 
suitable in this flow field because the assumption for the extra strain rate is no 
longer valid. In this case, the approximation procedure will be performed based 
on a two-dimensional streamline coordinate (s,n) system, where s-coordinate is 
in streamwise direction as i, j or k = 1 , and n-coordinate represents the direction 
normal to s-coordinate as i, j or k = 2. The velocity components are designated 
as U = Uj and u = u, in the streamwise direction, and V = U 2 and v = u 2 in the 
normal direction. The term— U/R is known to represent the extra rate of strain 
in the s-n coordinate system in addition to the main rate of strain 3U/3n (see 
Ref. [30]). By retaining the terms associated with these two strain rates and 
following the similar procedure as that in the Cartesian coordinate, the modified 
eddy viscosity fi t can be derived as shown in Appendix C, and hence 
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(3.22) 


k 2 2<() 

^‘ = P T- 


l - R f - 4> 


P r R f 2 + 4R f + 1 


1 -R, 


where 


R f = 


2 — 

R 


1 + 


n 

■r 




au + u 

an + r 


and R is the radius of surface curvature, n is the normal distance away from the 
wall. For convex curvature, the radius R has the same direction as n, so R is 


positive; while the radius R is in the opposite direction as n for concave 
curvature, hence R is negative. 

It is interesting to note that equation (3.22) is exactly the same as 
equation (3.19), except the definition of the flux Richardson number R f is 
different. This is very practical for engineering applications which implies no 
change of formulation for different coordinate systems. It is easy to recognize 
that the eddy viscosity in equation (3.22) will be diminished for convex 
curvatures due to the positive radius of surface curvature and flux Richardson 
number. The increase of the eddy viscosity of course will be introduced with 
the application of concave curvature because of the negative radius of surface 
curvature and negative R f . 
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3.3.3 Transverse curvature 


Transverse curvature occurs when flows have swirl. The cylindrical 
coordinate system, X ( = [x, y, 0], will be used with x = X, in longitudinal 
(axial) direction, r = X 2 in the radial direction, and 0 = X 3 in circumferential 
direction. The symbols U; = [U, V, W] and Uj = [u, v, w] represent the mean 
and fluctuating velocity components in the x, r and 0 direction respectively. It 
is to be noted that the additional terms introduced by the swirling component 
appear not only in the Reynolds-stress productions but also in the Reynolds- 
stress convections, which was proved by Johnson [84], and Humphrey and 
Chang [85]. The full Reynolds stress equations are complicated, as shown in 
Appendix D, hence some assumptions are necessary in order to simplify the 
modified model. Axisymmetric, thin-shear-layer, and weak swirl are therefore 
assumed for the proposing model, which lead to 3/30 * 0, 3/3r > 3/3x, and W/r 
< 3W/3r. The weak swirl approximation is also consistent with the constraint, 
isotropic turbulence, of most turbulence models including the k-e model. The 
anisotropic effect becomes more important as the swirl increases. With the 
above assumptions and following the similar derivation procedure of two 
proposed models, the modified eddy viscosity jx, can be derived as shown in 
Appendix D, and hence 
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(3.23) 




1 + 2R f 

t^r7 


where 



? W3W 
r dr 


f su T 

rawv 

l*J * 

[dr j 


(3.24) 


It is evident that the flux Richardson number R f is positive when the 
extra rate of strain, W8W/(r8r), is positive, and so the eddy viscosity p, is 
reduced by this effect, which can be seen from equation (3.23). Whereas, the 
enhancement of the eddy viscosity p t can be achieved by the introduction of the 
negative flux Richardson number, i.e. the extra rate of strain is negative. 
Comparing equation (3.23) with equations (3.19) and (3.22), we can easily find 
even though the eddy viscosity formulations are different, the effects of the 
extra rate of strains, embedded in the flux Richardson numbers, on the eddy 
viscosity |i t are essentially the same. 

The computational predictions by these three modified models will be 
compared with the experimental measurements, and will be further discussed 
in chapter V (Results and Discussions). 
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3.3.4 Wall function 

As indicated in the literature, [11], [13]-[15]; the near-wall turbulent 
structure is not significantly affected by the surface curvature. Hence, the 
standard wall function, suggested by Launder and Spalding [75] and used in 
TEACH-based program, will be applied in the present research. Detailed 
description is explained in Ref. [75] and Ref. [79]. 


46 


CHAPTER IV 


NUMERICAL METHOD 

4.1 Equation Integration 

A two-dimensional and a three dimensional Navier-Stokes flow solvers 
(FDNS-2D and FDNS-3D) embedded with the k-e model, developed by Y.S. 
Chen [80-82], are used to test the proposed model in the present study. The 
transport equations of the mean flow and turbulence model are transformed into 
the general curvilinear coordinates (^, q, O- The system of coupled 
transformed equations is discretized into a set of linearized algebraic equations. 
In the discretization process, several techniques are employed to stabilize the 
numerical integration and maintain the same order of accuracy. 

1) A time-centered (Crank-Nicholson) difference scheme is utilized 
for the temporal discretization. 

2) A second-order central difference method plus a damping term are 
used to manage the convection terms. The damping term could 
be either fourth order, deduced from second-order upwind 
differencing scheme, or second order, resulted from first-order 
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upwind differencing scheme. With the coefficient of the damping 
term being specified explicitly, the numerical stability of the 
matrix solver can be enhanced. 

3) The diffusion, production and other source terms are discretized 
with the second-order central differencing scheme to maintain the 
same numerical accuracy. 

The set of discretized linear equations is solved by an alternating 
direction linear-relaxation method (ADI) with a simplified predictor-corrector 
algorithm. In this pressure-based predictor-corrector solution procedure, an 
explicit fourth-order pressure smoothing term is added to the velocity-pressure 
coupled discrete equation, derived approximately from the discrete continuity 
and momentum equations, to inhibit the instability in the pressure solution. For 
each time step, the predictor-corrector loop will be executed iteratively until the 
conservation of mass is satisfied. 

4.2 Grid Generation 

Although there will be a variety of flow geometries, orthogonal grids 
system would be generated to avoid the ill-conditioned matrix caused by large 
grid skewness. Figure 1 illustrates the grid system for a backward-facing step, 
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and the mesh system for a 30°-45° S-bend duct is exhibited in Figure 2. Grid 
packing near the wall and large gradient zones may also be observed in Figures 
1 and 2. 
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CHAPTER V 


RESULTS AND DISCUSSIONS 


To evaluate the performance of the proposed turbulence model, various 
curved flows are examined in the present study. According to the previous 
classification of streamline curvature, numerical computations are performed on 
different geometries for each type of curvature effects. The numerical 
predictions and discussions are divided into 1) local curvature, 2) longitudinal 
curvature, and 3) transverse curvature. 

5.1 Local Curvature 

A backward-facing step with separation and reattachment provides the 
geometry for local flow curvature. Driver and Seegmiller [83] conducted an 
experimental investigation on the incompressible turbulent flow over a 
backward-facing step with an area ratio of 8:9. The grid system is shown in 
Figure 1 and the detailed infonnation about the test configuration is provided 
in Figure 3. To study the effect of varying modeling constant values (Ci and 
C 2 ) and the dependence of the cell size on the numerical predictions, two sets 
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of grid systems (61 x 41 or 1 1 1 x 45) along with different values for modeling 
constants (C, = 1.5 or C, = 1.8) are tested for the present model. 

First, a comparison is made between the present model and the standard 
k-e model with C, = 1.5 and 61 x 41 grids. The improvement by the present 
model is discernible from the locus of flow reversal illustrated in Figure 4. The 
result shows that the flow separation from the step wall generates the effect of 
convex curvature, which attenuates the eddy viscosity and actually causes the 
flow to reattach further downstream. Since there is no mechanism in the 
standard k-e model to simulate the curvature effect, a predicted convex shear 
layer exhibits higher viscosity and early reattachment takes place. Figure 5 
demonstrates better predictions by the present model on the streamwise velocity 
profiles. 

The more rapid change of the velocities in the recirculation zone 
predicted by the standard k-e model indicates the higher energy generated by 
the larger eddy viscosity. The turbulent kinetic energy profiles shown in Figure 
6, indicate the reduction of the turbulent kinetic energy by the convex curvature 
in the present model. However, both models underpredict the turbulent kinetic 
energy when the flow is near the reattachment location. This is possibly 
attributed to the effect of large-scale eddy, which becomes pronounced inside 
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the recirculation zone; and the current single-time-scale model fails to capture 
this phenomenon. 

The prediction on the Reynolds shear stress -puv by the present model, 
as shown in Figure 7, appears to be reasonably accurate. The success of the 
prediction of the Reynolds shear stress in the recirculation zone reveals validity 
of the eddy viscosity formulation. Meanwhile, it is worth noting that the 
present model costs 755 seconds of Cray X-MP CPU time (for 3122 time 
iterations) to converge to an accuracy of 5 x 10 4 , while 572 seconds of Cray 
X-MP CPU time (for 2486 time iterations) are required by the standard k-e 
model to converge to the same accuracy. It is encouraging that with all the 
improvements of the numerical predictions by the present model, only three 
minutes of Cray X-MP CPU time in addition are introduced. This shows the 
practicality and feasibility of the current approach for the complex engineering 
applications involving curvature. The velocity contours, static pressure 
contours, and velocity vector plots for both the standard k-e model and the 

present model are displayed in Figures 8-13. 

To demonstrate the insensitivity of the modeling constant C, value on 
computational results, a comparative investigation between C, = 1 .5 and C, = 
1.8 is performed. Figures 14-17 illustrate almost identical results of the flow 
reversal locations, streamwise velocity profiles, turbulent kinetic energy profiles, 
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and Reynolds shear stress profiles between the two C, values. The above 
exercise satisfies our previous statement that the performance of the present 
model is insensitive to a reasonable variation of C, value. The remaining 
computations are performed with the empirical coefficient C, specified as 1.5. 

By increasing the grid numbers from 61 x 41 to 111 x 45, a better 
description of not only the change of streamline curvature but also the discrete 
flow field is expected to be achieved. From Figures 18-21, it appears that the 
numerical calculations of the locus of flow reversal, streamwise velocity, 
turbulent kinetic energy, and the Reynolds shear stress profiles are independent 
of the grid size. With the same accurate results, the medium grids instead of 
finer grids will be employed to save the computing time. 

To achieve further verification of the proposed model with local 
curvature effect, the numerical prediction on a backward-facing step with an 
area ratio of 2:3, tested by Kim, Kline, and Johnston [86], is conducted and 
the computational result is compared with the experimental data. The layout 
of the backward-facing step is sketched as shown in Figure 22 and will be 
discretized into an 85 x 51 mesh system. The improvement of the present 
model, as indicated in Figure 23, is consistent with the result of the previous 
backward-facing step case. It is evident that the locus of flow reversal, plotted 
in Figure 23, reveals the superiority of the proposed model over the standard 
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k-e model in predicting the larger recirculation zone. This is indicative of the 
attenuation of the eddy viscosity caused by the effect of convex curvature. The 
measured reattachment length x R is 7±1 step heights; and it verifies the 
prediction of the present model with later reattachment. The variation of x R is 
due to the inherent unsteadiness of the flow field. 

Figure 24 exhibits better performance of the present model in capturing 
the streamwise velocity profiles. The reduction of the turbulent kinetic energy 
by the effect of convex curvature in the recirculation zone is substantiated by 
the agreement between experimental data and numerical predictions of the 
present model, which is illustrated in Figure 25. The Reynolds shear stress -uv 
profiles shown in Figure 26 reveal better predictions of the present model. The 
successful predictions of the turbulent kinetic energy and the Reynolds shear 
stress by the present model indicates that the proposed formulation for the eddy 
viscosity is capable of describing the characteristics of the flow separation with 
curvature effect. Figure 27 also illustrates that the proposed model can predict 
the wall static pressure coefficients better than the standard k-e model. 

Although the overall predictions of the present model are reasonably 
successful, there is some discrepancy between experimental results and 
numerical predictions. The disagreement can be attributed not only to the 
unsteadiness of the flow field but also to additional turbulence phenomena, such 
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as large-scale eddy and anisotropy not accounted in our theory. Therefore, it 
is suggested that further improvement on the proposed model be made by taking 
into account the physics of the large-scale eddy and anisotropy. This will be 
explained further in chapter VI-NUMMARY" of this report. The detailed 
descriptions of the flow field computed by both the standard k-e and the present 
model, such as velocity contours, velocity vectors, and particle traces, are 
shown in Figures 28-31. 

5.2 Longitudinal Curvature 

The most typical longitudinal curvature occurs in a flow through a 
curved duct. An S-bend duct is a more complicated geometry because the 
second bend of the duct creates an effect opposite to that of the first bend. In 
this case, if the model can not capture well the effect of curvature in the first- 
bend section, the error will accumulate in the second bend and cause the 
collapse of the model. A 30°-45° S-bend duct with an aspect ratio of 5.6, 
measured by Butz [23], will be tested against the present model. The detailed 
description of the facility and the inlet flow conditions are shown in Figure 32. 
An 88 x 51 grid system is constructed, as shown in Figure 2; and the inlet 
velocity and turbulent intensity profiles are prescribed. 


55 



The comparison of static pressure coefficients along the curved wall is 
illustrated in Figure 33. It seems that in the first bend, both models fail on the 
upper surface (concave side), but succeed on the lower surface (convex side). 
In the second bend, the present model does a better job in predicting surface 
pressure than the standard k-e model. One might question the success of the 
predicted results, on the upper surface (convex side) in the second bend by the 
present model. However, as the measured data exhibits some randomness and 
the present model predicts the exit static pressure relatively well, the numerical 
prediction of the present model can be considered to be reasonably successful. 
The discrepancy on the upper surface in the first bend is suspected of being 
caused by the Gortler-like secondary flows. This occurs at the onset of the 
concave curvature, where the normal velocity is numerically set equal to zero 
numerically at the entrance plane. There is no experimental measurement of the 
normal velocity component at the entrance plane. The error of predicted static 
pressure on the upper surface (concave side) in the first bend is considered to 
be the primary source of deviation of the numerical calculation from 
experimental result on the same surface (convex side) in the second bend. 

The longitudinal velocity profiles at the inlet, inflection and exit planes 
from the experimental results and both computational models are compared in 
Figures 34-36. The results demonstrate that the proposed turbulence model has 
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a better agreement with the experimental data than the standard k-e model. The 
slight deviations of the longitudinal velocity from the measured data at the 
inflection and exit planes are attributed to the effect of the small secondary 
flow. With 88 x 51 grids, the present model converges to 5 x 10' 4 in 110 
seconds of Cray X-MP CPU time with 199 time iterations, while the standard 
k-e model consumes 116 seconds with 218 time iterations. It is interesting to 
note that the present model has improved the numerical predictions without 
increasing the computing time substantially. The contour plots of x-component 
velocity, y-component velocity, and static pressure are also shown in Figures 
37-39 for both the present and the standard k-e models. 

The 30°-45° S-bend duct does not demonstrate the superiority of the 
proposed model over the standard k-e model very well because the surface 
curvature is relatively small. The centerline radius of curvature is five times the 
duct width, and hence the effect of streamline curvature does not dominate the 
flow field significantly. A 2-D curved duct with small radius of curvature 
(strong curvature), therefore, is required to evaluate the performance of the 
present model in predicting the development of turbulence structures subjected 
to strong curvature effect. A 10-by-100 centimeter, 1 80-degree-tum water 
tunnel, with a 10 cm centerline radius of curvature, investigated by Sandborn 
and Shin [87], is employed as the next test case. The configuration of the 180 5 
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turn-around duct (TAD) and inlet conditions are described in Figure 40. The 
flow field in the 180° turn-around duct with strong curvature, especially, is of 
great engineering interest due to its resemblance to the flow passage 
downstream of the turbine in the fuel prebumer of the Space Shuttle Main 
Engine (SSME) as well as STOVL applications. Complex shear layers, regions 
of separation, high levels of unsteadiness, and three-dimensional structure may 
occur in this type of flow field. A 141 x 41 orthogonal grid system is 

constructed for this case, as shown in Figure 41. 

The longitudinal velocity profiles in Figure 42 display the similarity 
between the predictions of the present model and those of the standard k-c 
model upstream of the flow separation. However, the results demonstrate the 
success of the present model and the collapse of the standard k-e model in 
capturing the separation bubble on the inner surface (convex curvature) near the 
180° turn. Once again, as with its failure to predict the flow in a backward- 
facing step, the standard k-e model is inadequate in describing the change of 
the eddy viscosity (or the length scale) attributed to the effect of streamline 
curvature. The larger eddy viscosity along the inner surface predicted by the 
standard k-e model prevents the occurrence of flow separation; whereas, smaller 
eddy viscosity computed by the present model, benefiting from curvature 
correction, enables the natural onset of flow separation. The numerical results 
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of the present model reveal that flow starts to separate after 170° turn. Figure 
43 provides a detailed description about the longitudinal velocity distribution in 
the separation bubble downstream of 180° turn. The results indicate that the 
proposed implementation to the standard k-e model is necessary and proper. 

The wall static pressure coefficient C p is plotted as shown in Figure 44. 
It is surprising that the standard k-e model predicts much better than the present 
model even without predicting flow separation. However, the pressure 
distribution calculated by the present model does show the signature of flow 
separation where the static pressure recovery along the inner surface is retarded. 
That is the location where the predicted static pressure coefficient C p by the 
present model starts to deviate from the measured data. Hence, the relative 
value of the measured static pressure is doubtful. It was later discovered in 
Ref. [87] that screens were placed at the outlet of the channel, which is about 
four channel widths downstream of 180° turn, to increase the flow exit 
resistance such that the water completely filled the channel. However, the 
screen raised the adverse pressure gradient at the exit and caused the exit static 
pressure to be larger than it would have been. To support this argument, 
another 180° turn-around duct will be examined next. Some information about 
the flow field such as x- and y-component velocity contours, static pressure 
contour, velocity vectors, and particle traces are shown in Figures 45-49. The 
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separation bubble can be observed clearly from the plots of the velocity vectors 
and particle traces. 

A further validation will be made by examining a 2-D U-duct 
investigated by Monson and Seegmiller [88]. The geometry of the coordinate 
system, and the inlet conditions, are illustrated in Figure 50. The ratio of the 
centerline radius of curvature to the channel width is equal to unity, while a 234 
x 101 mesh system is constructed due to longer inlet and exit ducts as well as 
higher Reynolds number. 

The numerical results reveal that both the present and standard k-e 
models predict flow separations. Nevertheless, the flow field calculated by the 
present model is observed to separate earlier (ahead of 1 80° turn) and reattach 
later, which is consistent with the experimental measurement. A later flow 
separation (preceded by 180° turn) and an earlier reattachment are suggested by 
the standard k-e model. This is confirmed by the comparison of longitudinal 
velocity profiles as shown in Figure 5 1 . Nevertheless, the present model seems 
to under-estimate the thickness of the separation bubble. It is then found in 
Ref. [88] that the thickness of the separation bubble in the Re = 10 6 case is 
larger than that in the Re = 10 5 case; while most of the turbulence models, 
including the present model, predict it differently. From the physical 
arguments, the Reynolds number trend in experimental results seem to be doubtful. 
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The distribution of wall static pressure coefficient C p , plotted in Figure 
52, supports the previous argument. The present model has better agreement 
than the standard k-e model in predicting the static pressure. Furthermore, the 
signature of flow separation presented by the proposed model does match the 
flow characteristic near the exit of 180° turn. However, both models fail to 
predict the static pressure at the downstream exit plane. It seems that the actual 
flow separates earlier than the present model, and this causes the present model 
to predict larger static pressure recovery at the exit. It is apparent from these 
comparisons that the present model performs a better overall prediction, due to 
the adjustment to the effect of surface curvature, than the standard k-e model, 
especially on the convex surfaces where the attenuation of the eddy viscosity 
plays an essential role in the occurrence of flow separation. The contour plots 
of x- and y-component velocities as well as static pressure are illustrated in 
Figures 53-55, respectively. The plots of velocity vectors and particle traces, 
as shown in Figures 56-57, indicate an earlier, larger flow separation by the 
present model and a later, smaller flow separation by the standard k-e model. 

5.3 Transverse Curvature 

One representative of the class of flows with transverse curvature is the 
swirling flow. There are two constraints in selecting the test cases to validate 
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the proposed model with the effect of transverse curvature. First, the flow must 
have small swirl number such that the characteristics of the flow is consistent 
with the assumptions, i.e. weak swirl and isotropic turbulence, made in the 
present model. Second, confined swirling flows in straight pipes would be 
employed for the purpose of avoiding the interference from other curvature 
effects such as surface curvature effect, and also maintaining the axisymmetry 
of the flow. The numerical computations of the standard k-e and the present 
models for the effect of transverse curvature will be conducted in the form of 
2-D axisymmetric instead of truly three-dimensional test case. 

An experimental study on a confined swirling coaxial jet in a straight 
pipe with the swirl number equal to 0.21, investigated by Roback and Johnson 
[89], is therefore chosen to be the test case. A 151 x 45 grid system is built up 
for this flow geometry, and the schematic of the configuration along with the 
inlet conditions are illustrated in Figure 58. The inlet velocity profiles and 
turbulence quantities are specified at 5 mm downstream of the jet nozzle. With 
such a small swirl number, a vortex breakdown occurs in this flow field 
because of the combination of axial and radial pressure gradients. The radial 
pressure gradient is induced by the centrifugal force, while the axial pressure 
gradient is attributed to both centrifugal force of the swirling flow and the 
diffusing geometry of the experimental set up, i.e. dump diffuser. 
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Consequently, the comparison of the mean axial velocity along the centerline, 
demonstrated in Figure 59, is one indication for the performance of both 
models. Although both models fail to depict the size of the central recirculation 
zone, as caused by overpredicting the eddy viscosity, the present model appears 
as an improvement over the standard k-e model. The large eddy viscosity 
computed by both models is suspected to be strongly related to the specification 
of the length scale at the inlet boundary. The inlet length scale is not available 
from the experimental results, and is assumed to be one percent of the pipe 
diameter. The change of the inlet length scale is expected to greatly alter the 
numerical results; but it is not the main concern in the present study. 

A detailed comparison of the mean axial velocities is illustrated in Figure 
60. The experimental results show some random scatter, which is partially 
caused by the non-axisymmetric phenomenon. It is evident that the present 
model does provide minor improvements over those by the standard k-e model, 
but the improvements are insufficient to satisfy the physical characteristics of 
the flow. In addition to the uncertainty of the inlet length scale, the following 
reasons may be cited: 1) the additional streamline curvature caused by the 
central recirculation zone affects the flow field, 2) asymmetric flow, 3) 
anisotropic turbulence characteristics, and finally, the inadequacy of the current 
implementation approach for swirling flows. To identify the causes for only a 
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minor improvement by the present model, another confined swirling flow case 
will be examined later on. 

The mean azimuthal and radial velocity profiles, plotted in Figures 61-62, 
indicate that the present model predicts better than the standard k-e model, but 
the results are relatively similar. A detailed description of the flow field such 
as the axial and the radial velocity contours, the velocity vectors, and the 
particle traces are reported in Figures 63-66. The central recirculation zone 
predicted by the present model is larger than that suggested by the standard k-e 
model, and can be clearly observed from the particle trace plot. 

The next test case for the transverse curvature effect is the confined 
swirling flow studied by Weske and Sturov [90]. The layout of the test 
section and the specification of the inlet conditions are exhibited in Figure 67 . 
With the inlet velocity profiles and turbulence quantities prescribed at 3.5 cm 
downstream of the swirl generator, the inlet swirl number is calculated to be 
0.43 and has a distribution of solid body rotation. 

The axial velocity profiles are compared in Figure 68 but with the 
absence of experimental data, because it is not reported in Ref [90]. As shown 
in Figure 68, almost identical numerical results are predicted by both models. 
In Figure 69, the comparisons of the circumferential velocities illustrate that the 
standard k-e model predicts slightly better than the present model. This is 
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contrary to the results in the previous swirling flow case. However, both 
models predict faster mean flow decay of the swirling flow than the 
experimental results indicate. The present implementation approach to include 
the effect of transverse curvature (swirling effect) is proved to be inadequate 
based on the results of two swirling flows. It is believed that the derivation of 
the eddy viscosity ji, from the equation of the Reynolds shear stress uv is 
improper. The reason for the deficiency is that the other Reynolds shear stress 
components, uw and vw, are of the same order of magnitude as the shear stress 
tensor uv in the swirling flows. Therefore, difficulty will be encountered in 
choosing the Reynolds shear stress component from which the scalar eddy 
viscosity would be determined. Employment of the scalar eddy viscosity is not 
considered to be suitable for the flow with swirling effect. Instead, a vector 
eddy viscosity for each Reynolds shear stress, or calculating each Reynolds 
stress tensor through algebraic equations, as derived similarly in the present 
study, is recommended for this type of flow for future studies. 
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CHAPTER VI 


SUMMARY 


6.1 Conclusions 

The results of the test cases, employed to study the effects of local and 
longitudinal curvatures, show the success of the proposed eddy viscosity 
formulation in improving the numerical prediction capability of the standard k-e 
model consistently. Moreover, the present model offers the following 
advantages, namely 1) generality— same formulation for the eddy viscosity in 
both Cartesian and Streamline coordinates, 2) convenience— the formulation for 
the eddy viscosity is easy to be adopted into various one- and two-equation 
models, 3) efficiency— the required computer CPU time by the present model 
is nearly the same as the standard k-e model. 

In the backward-facing step flows, the present model is indeed capable 
of dictating the effect of local curvature rendered by flow separation and hence 
capturing the attenuation of the eddy viscosity. The present model, however, 
still underpredicts the reattachment point. This problem is attributable to the 
effects of large-scale eddy and anisotropy in the recirculation zone, especially 
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the latter one. The flow inside the recirculation zone is basically dominated by 
a large-scale eddy with low Reynolds number. This is a deficiency of the high- 
Reynolds-number k-e model used in this study. It is also known that the 
anisotropic effect becomes dominant as the flow approaches the wall, thus the 
effect plays an important role near the reattachment region. Consequently, the 
proposed model could be significantly improved with the consideration of these 
two effects. 

In the curved-duct flows, the present model demonstrates a good 
agreement with the measured data. The major success of the present model is 
to predict the flow separation in curved ducts by accounting for the effect of 
wall curvature. With the absence of the curvature effect, the standard k-e 
model appears to predict the flow without separation or with later and smaller 
separation zone. However, the discrepancy between the measured data and the 
numerical results takes place when the present model estimates the size of 
separation bubbles. As discussed in chapter five, the experimental results are 
believed to be questionable. Hence, a more detailed and reliable experiment for 
the 180 -degree turn- around duct is required to further validate the present model 
in predicting the size and location of the separation bubble. 

For the flows with swirling effect, the present model does not improve 
the numerical predictions consistently. In spite of depicting the size of vortex 
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breakdown in the confined swirling coaxial jets better than the standard k-e 
model, the present model overpredicts the eddy viscosity and results in faster 
decay of the circumferential velocity in the confined swirling pipe flow. The 
current approach to include the effect of transverse curvature, deducing the 
scalar eddy viscosity from the equation of the Reynolds shear stress uv, is thus 
considered to be improper for the flows subjected to swirling effect. With the 
dominance of the anisotropic effect in the swirling flows, each Reynolds shear 
stress term can be of the same order of magnitude and should have its 
corresponding eddy viscosity formulation, respectively. Therefore, the effect 
of transverse curvature should be accounted for either by employing a vector 
form of the eddy viscosity or by solving for each Reynolds stresses directly 
through approximated algebraic equations. Although, the numerical predictions 
of the swirling flows by the present model show deficiencies, a qualitative 
improvement is accomplished, which reveals that the extra rate of strain should 
be included to account for the effect of transverse curvature. The 
implementation of transverse curvature is expected to be more complicated than 
those for the local and longitudinal curvatures due to turbulence anisotropy. 

The present model demonstrates improvements to the prediction 
capability of turbulent flows dominated by local and longitudinal curvatures 
over the standard k-e model. The effect of transverse curvature, as in swirling 
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flows, requires more elaborate formulation of the eddy viscosity and should be 

included in future work. 

6.2 Recommendations for Future Work 

1) The present model demonstrates good predictions in the backward-facing 
step cases, where the effect of local curvature prevails. It is 
recommended that the anisotropic effect be included, hence the 
reattachment length can be predicted more precisely. 

2) A more extensive experimentation on the 180-degree tum-around duct 
is required to evaluate the performance of the present model in 
predicting the separation bubble rendered by the effect of wall curvature 
(longitudinal curvature). 

3) A vector form of the eddy viscosity or solving the algebraic Reynolds 
stress equations should be employed to account for the effect of 
transverse curvature (swirling flows). The effect of strong anisotropy 
will result in the breakdown of the algebraic eddy viscosity approach. 
The adoption of the extra rate of strain induced by the curvature effect 
should still be valid, but each components of the vector eddy viscosity 
should be attained from corresponding Reynolds stresses. 
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Figure 1 Grid System for a Two-Dimensional Backward-Facing Step Geometry 



Figure 2 Mesh System for a Two-Dimensional S-Bend Duct Geometry 
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X = 0 


Tunnel Geometry: h = 1.27 cm 

W in = 8h 


L = 30h 
L u = 4h 

aspect ratio = 12 


Inlet Conditions: U ref — 44.2 m/sec 

0.128 
8b L = 1.9 cm 
Re e = 5000 


Figure 3 Backward-Facing Step Flow Experimental Geometry and Inlet 
Conditions by Driver and Lee Seegmiller [83] 
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Data : Driver and Lee Seegmiller [ 83 ], : Standard, : Present (C 
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(61 x 41 grids) 


Data : Driver and Lee Seegmiller t 83 ], : Standard, : Present (C 



Continued 



Data : Driver and 
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Figure 6 Turbulent Kinetic Energy Profiles in a Backward-Facing Step Flow (6 
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Data : Driver and Lee Seegmiller [ 83 ], : Standard, : Present (C 



87 


Figure 7 Continued 


(a) present model 



Figure 8 


Streamwise Velocity Contour in a Backward-Facing Step Flow 
(61 x 41 grids) 



Figure 9 



(b) standard k-e model 


Streamwise Velocity Contour near the Backward-Facing Step 
(close view) 
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(b) standard k-e model 


Figure 10 Transverse Velocity Contour in a Backward-Facing Step Flow (61 
x 41 grids) 
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(a) present model 




Figure 11 Transverse Velocity Contour near the Backward-Facing Step 
(close view) 
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(b) standard k-e model 


Figure 12 


Static Pressure Contour in a Backward-Facing Step Flow (61 x 41 
grids) 
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(a) present model 



(b) standard k-e model 


Figure 13 Velocity Vectors in a Backward-Facing Step Flow (61 x 41 grids) 
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Figure 14 Locus of Flow Reversal in a Backward-Facing Step Flow with 
Different Modeling Constants 



Data : Driver and Lee Seegmiller [ 
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Figure 15 



Data 
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Continued 



Bad 



Data 





Data : Driver and Lee Seegmiller [ 83 ], 



Figure 17 Continued 



Figure 18 Locus of Flow Reversal in a Backward-Facing Step Flow with 
Different Grid Size 
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Data : Driver and Lee Seegmiller [ 83 3, : 101 



Figure 19 Continued 



Data : Driver and Lee Seegmiller t 83 J, 
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Figure 20 Turbulent Kinetic Energy Profiles in a Backward-Facing Step Flow with Different Grid Size 


Data : Driver and Lee Seegmiller t 83 ], : 101 
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Data : Driver and Lee Seegmiller [ 83 ], 
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Figure 21 Continued 
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Figure 22 Backward-Facing Step Flow Configuration and Inlet Conditions [86] 



Figure 23 Locus of Flow Reversal in a Backward-Facing Flow (85 x 51 
grids) 
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Data : Kim, 



Figure 24 





Data : Kim, 
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Data : Kim, Kline, and Johnston [ 86 ], : Standard, : Present (C 
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Data : Kim, Kline, and Johnston [ 86 ], : Standard, : Present (C 
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Data : Kim, Kli 
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Figure 26 




Data : Kim, 
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Figure 26 
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x / h 


Figure 27 Static Pressure Coefficient Distributions along the Wall in 
Backward-Facing Step Flow (85 x 51 grids) 
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Figure 30 Velocity Vectors in a Backward-Facing Step Flow [861 (85 x 51 grids) 





(a) present model 



(b) standard k-e model 

Figure 31 Particle Traces in a Backward-Facing Step Flow (close view) 
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Re w = 131200 
U, ve = 63 ft/sec 
Aspect ratio = 5.6 

Figure 32 A 30°-45° S-Bend Duct Geometry and Inlet Conditions by L.A. 
Butz [23] 


Pressure Coefficient C 



Figure 33 Static Pressure Coefficient Distribution for a 30°-45° S-Bend Duct 
(88 x 51 grids) 
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Figure 34 Longitudinal Velocity Profile at the Entrance Plane of a 30°-45° 
S-Bend Duct (88 x 51 grids) 
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Figure 35 Longitudinal Velocity Profile at the Inflection Plane of a 30°-45° 
S-Bend Duct (88 x 51 grids) 
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Figure 36 Longitudinal Velocity Profile at the Exit Plane of a 30°-45° S- 
Bend Duct (88 x 51 grids) 
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Figure 37 



(b) standard k-e model 


X-Component Velocity Contour in a 30°-45° S-Duct (88 x 51 
grids) 
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(a) present model 



(b) standard k-e model 


Figure 38 Y-Component Velocity Contour in a 30°-45° S-Duct (88x51 grids) 
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(a) present model 



(b) standard k-e model 


Figure 39 Static Pressure Contour in a 30°-45° S-Duct (88 x 51 grids) 
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Test Section Geometry Inlet Conditions 


W = 10 cm 
R ; = 5 cm 
R c = 10 cm 
R 0 = 15 cm 
L u = 17.27 cm 
L d = 43.2 cm 
aspect ratio =10 


U ref =281.9 cm/sec 
Re w = 213000 


Figure 40 The Geometrical Sketch and Inlet Conditions for a 180° Turn 
Around Duct by Sandbom and Shin [87 J 



Figure 41 


Grid System for a Two-Dimensional 180° Turn-Around Duct 
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Figure 44 Wall Static Pressure Coefficient Distribution in a 180° Turn- 
Around Duct [87] 
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Figure 45 X-Component Velocity Contour in a 180° Turn-Around Duct [87] (141 x 41 grids) 
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Figure 46 Y-Component Velocity Contour in a 180° Turn-Around Duct [87] (141 x 41 grids) 
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Figure 47 Static Pressure Contour in a 180° Turn-Around Duct [87] (141 x 41 grids) 



e 48 Velocity Vectors in a 180° Turn-Around Duct [87] (141 x 41 grids) 
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Figure 49 Particle Traces in a 180° Turn-Around Duct [87] (141 x 41 grids) 



Figure 50 The Configuration and Inlet Conditions for a 2-D U-Duct [88] 



Figure 





Pressure Coefficient C 



Figure 52 Wall Static Pressure Coefficient Distribution in a U-Duct (234 x 
101 grids) 
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(b) standard k-e model 

Figure 53 X-Component Velocity Contour in a 2-D U-Duct [88] 
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(b) standard k-e model 

Figure 54 Y-Component Velocity Contour in a 2-D U-Duct [88] 
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(b) standard k-e model 

Figure 55 Static Pressure Contour in a 2-D U-Duct [88] 
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(b) standard k-e model 
Figure 57 Particle Traces in a 2-D U-Duct [88] 
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Illustrations of Test Section, Coordinate System, and Inlet Conditions for a Confined Swirling Jet [89] 



(m/sec) 



Figure 59 Mean Axial Velocity along the Centerline in a Confined Swirling 
Jet (151 x 45 grids) 
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Data : Roback and 



Figure 61 







k-e model 


in a Confined Swirling Jet [89] 
56 



(b) standard k-e model 


Figure 64 Radial Velocity Contours in a Confined Swirling Jet [89] 
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(a) present model 



(b) standard k-e model 


Figure 65 Velocity Vectors in a Confined Swirling Jet [89] 
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(a) present model 



Figure 66 Particles Traces in a Confined Swirling Jet [89] 
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Figure 67 Schematics of the Test Section and the Inlet Conditions for a Confined Swirling Flow [90] 


Data : Weske and Sturov [ 90 ], : Standard, : Present (C 



Figure 68 Axial Velocity Profiles in a Confined Swirling Flow (301 x 41 grids) 
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APPENDIX B 


DERIVATION OF NEW FORMULATION FOR EDDY VISCOSITY IN 

CARTESIAN COORDINATES 


The algebraic equation for the Reynolds stress u ; Uj in equation (3.17) is 


k e 


P - Is. P 
o a u r 



where P^ is the production rate of the Reynolds stress given by 


(B.l) 



au. au 

U.U._— i + u : u. 

* k ax. 




ax. 


(B.2) 


and P r is the production rate of the turbulent kinetic energy given by 



au* 

ax k 


(B.3) 


For two-dimensional Cartesian (x,y) coordinates, the streamline direction 


(x) will be in accordance with the indices i, j and k equal to 1, while the indices 
i, j and k=2 represent the transverse direction (y). Therefore, the above 
equations can be written as 
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3x 3y 3x 3y 
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Following Bradshaw’s observation [30]-[31], only the terms associated with 


3U/3y and 3V/3x will be kept, and hence 
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dy 
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au av 
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(B.14) 


According to Bradshaw’s analogy, the flux Richardson number R f is the 
negative ratio of the turbulent kinetic energy production in the transverse 
direction to that in the streamwise direction, i.e. 

_ - p 22 _ ( -2uv3V/3x ) = -3V/3x (B.15) 

f " ~P^~ (-2U73U/3y) ^U/3y 

Combining equation (B.15) with equation (B.14), it yields 

= p _i_ ( B16 ) 

3y r 1 - R f 

= -P -A- ( B - 17 > 

3x r 1 - R f 

Substituting equation (B.16) into equation (B.ll) and equation (B.17) into 
equation (B.12) will give 

From equations (B.4) and (B.18) as well as equations (B.5) and (B.19), it can 
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be shown that 
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Equation (B.13) also can be rewritten as 
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ay 
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A new form for the equation of P 12 can be obtained with the help of equation 


(B.20) and (B.21) 
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According to equation (B.6), the algebraic equation for the Reynolds shear 
stress can then be obtained by incorporating with equation (B.23) as 


B4 


-uv = 


k 2 au 2 

e dy 3 



R f 2 + 4R f + 1 
»-Rf 


(B.24) 


B5 


APPENDIX C 


DERIVATION OF NEW FORMULATION FOR EDDY VISCOSITY IN 

STREAMLINE COORDINATES 


Following Rodi’s [61] step, the algebraic equation for the Reynolds stress 


UjUj is obtained as 
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where Pjj, the production rate of the Reynolds stress is 
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and P r , the production rate of the turbulent kinetic energy is 


au. 


P = -u.u. 

r 1 k ax. 


(C.3) 


With the selection of a two-dimensional streamline (s,n) coordinates, the 
streamwise direction (s) will be aligned with the indices i, j and k equal 1, 
while the indices i, j and k equal to 2 represent the normal direction (n). 
Therefore, the above equations can be expressed as 
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where R is the radius of surface curvature and n is the normal distance away 
from the wall. 

Since the dominant extra rate of strain is U/R, only those term coupled 
with U/R and the main strain au/3n will be remained, and hence 
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According to Bradshaw’s analogy [30]-[31], the flux Richardson number R f is 
the negative ratio of turbulent energy in the normal direction to that in the shear 
direction, i.e. 
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Combining equations (C.14) and (C.15) will yield 
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The above equations show that equations (C.ll) and (C.12) can be expressed 


as 
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P n = 2P 
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(C.18) 
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Substituting equation (C.18) into equation (C.4) and equation (C.19) into 
equation (C.5) will obtain 
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Since the equation (C.13) can be rewritten as 
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a new expression for P 12 can be deduced with the aid of equations (C.20) and 
(C.21), which is 
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By putting equation (C.23) into equation (C.6), an algebraic equation for the 
Reynolds shear stress can be resolved as 
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The Boussinesq’s eddy-viscosity reveals that 


-puv = ji t 


/i + Oau u 


V 


R 


dn R 


so that we can obtain eddy viscosity p t as 


(C.24) 


(C.25) 


k 2 2<}> 

n, = p — -r 

e 3 


P R f 2 + 4R f + 1 

1 - R - 4, _L _[ 

r E 1 - R, 


(C.26) 
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APPENDIX D 


DERIVATION OF NEW FORMULATION FOR EDDY VISCOSITY 
WITH TRANSVERSE CURVATURE EFFECT 


With the adoption of the cylindrical coordinates, the coordinate axis, 
mean and fluctuating velocity components can be expressed in the tensor 
notation as X ; = [x, r, 0], Uj = [U, V, W], and u ; = [u, v, w]. According to Ref. 
[84] and Ref. [85], the transport equations of Reynolds stresses UjUj and 
turbulent kinetic energy k can be expressed as follows: 


U. = -2uu. + Ps + D - e n 

k ax. k ax. 11 11 " 


(D.l) 


a^ 7 _ w— 

U. _ — - 2 — vw = -2 

ax k r 


r 


vu„ 


av w 


ax k r 


. VW 


+ ^ S 22 + ^22 e 22 


(D.2) 


TT dw 1 0 W — ~ 

U. + 2 — vw = -2 

k ax. r 


wu,. 


aw v 


ax k r 


. w 


+ Ps 33 + D 33 e 33 


(D.3) 
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a uv w — 

U. - — - — uw 
k ax k r 


— au — av 

vu, - — + uu. 


ax. 


ax. 


w — 

uw 

r 


T , a U w w — 

U. _ — + — uv 

k ax, r 


+ Ps ]2 + D ]2 e ]2 


— au — aw 

wu. _ + UU. _ 

k ax k k ax K 


+ Ps, 3 + D ,3 e, 3 


V — 
— uw 

r 


T T avw w / -j — 7 - 
U, - — + — ( v - w z ) = - 


ax k r 


(D.4) 
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av — aw 

+ vu k 

ax, k ax. 


— w 

. vw w 


(D.6) 


+ PS 23 + D 23 - e 23 


U, 


ak 


k ax. 


( — au w — v ^ 

U U. - VW + — w 

'* ax k r r , 


+ P k + D k " 6 


(D.7) 


where subscripts i and k are the tensor indices, ax k = [3x, 3r, ra0], and Ps iJS D^, 
and are the pressure-strain, diffusion, and dissipation of Reynolds stresses, 
respectively; while P k , D k , and e are the pressure-strain, diffusion, and 
dissipation of the turbulent kinetic energy. 

Therefore the production terms of Reynolds stresses and turbulent kinetic 
energy can be written as 


D2 



P = P = -2 

1 22 r rr ^ 


UV 


3 V — 2 3 V — 
+ V + vw 


3x 


3r 


( av _ 9 w ^ 

r50 r 
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P = P = -2 

33 r 09 * 
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UW + VW +• W + W + vw 

5x ar ra0 r r 
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12 Xf 
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- UV _ — - \ - VW - u 
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au 

ar 
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1 13 1 x6 
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( av + aw _ w "i 

r30 3r r 
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r30 dx 

V 

with the assumption of weak swirl, axisymmetric and slender shear layer, which 
means W/r < aw/3r, d/dQ ~ 0, and d/dr > a/3x, equations (D.8)-(D.14) can be 
simplified as 


D3 


(D.15) 


D D -> — dU 

P,i = P„ = -2uv 

u X* 5r 


W 
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r 23 r rQ V 


aw 
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The flux Richardson number R f of three-dimensional flows i 
Bradshaw [30] as 


£ extra V 2 production 

sum of u 7 and w 7 production 
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(D.16) 

(D.17) 

(D.18) 

(D.19) 

(D.20) 

(D.21) 
defined by 

(D.22) 



From equations (D.21) and (D.22), the production of turbulent kinetic energy 
can be derived as 


W 1 - R f 

P = -2vw 

r r R f 


(D.23) 


and so a new expression for can be obtained as 


R f 

P„ = -2P — 

22 r 1 - R, 


(D.24) 


According to equation (3.17), which is 


W = i 

k e 


P . - _8..P 

.j 3 ,j r 


lb-. 

3 ,J 


(D.25) 


an expression for the shear stresses uv and v 2 can be denoted as 


k i 


12 


(D.26) 
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v" _ <j) 


P 22 - 1 P 
22 3 r 


(D.27) 


Substituting equations (D.23) and (D.24) into equation (D.27) will yield 


— j i 

v z = _ k 
3 


P 1 + 2R f 

1 - <t> — L ! 

c 1 - R f 


(D.28) 


With the aid of equation (D.28), equation (D.18) can be rewritten as 
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(D.29) 


P 12 = ~-k 

12 3 




P 1 +2R, 


e 1 -R, 


au 

~aT 


Therefore, the Reynolds shear stress uv can be obtained from equations (D.26) 
and (D.29), which is 


— k 2 
uv = - — 

e 3 


2 <}> 


i m Pf 1 + ^ Rf ^ 
1-4) — — 

e 1 - R f 3r 


and also the eddy viscosity could be shown as 




1 + 2R f 

t ~ r ; 


(D.30) 


(D.31) 


The final step is to further simplify the expression of the flux Richardson 
number R f . Due to the assumption of weak swirl, i.e. W/r < 3W/3r, equation 
(D.22) becomes 



f _3U — 3W 

UV + VW 

3r 3r 

From equation (D.25), it can be seen that 


(D.32) 


uv _ P i 2 (D.33) 

VW P 2 3 
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Substituting equations (D.18) and (D.20) into the above equation will evolve the 
following expression 


uv _ dU/dr (D.34) 

^ 3W /dr 

Hence, the flux Richardson number R f can be derived from equations 
(D.32) and (D.34), which is 

2 w aw 

(D.35) 
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